A previous study had conducted both gene expression and metabolomics profiling of tissue samples from breast cancer patients. This vigenette will highlight the analysis we conduct on the breast cancer data.
IntLIM, available from Github, can be installed as in the documentation. Once IntLIM is installed, what is necessary is loading in the library.
rm(list = ls())
library(IntLim)
For breast cancer study, both gene expression (available on Gene Expression Omnibus Accession Number: GSE27751) and metabolomics (http://www.jci.org/articles/view/71180/sd/2) data are available online. Much of this data had been processed as previously described (3). Probes from the Affymetrix data not mapping to a gene symbol were removed. Additionally, as with NCI-60 data, only the probe corresponding to the highest mean expression was used for analysis when multiple probes corresponded to a single gene. This resulted in a total of 20,254 genes for 108 patient samples. The Metabolon data did not need to be filtered by coefficient of variation, as there were no technical replicates. The resulting data consisted of 536 metabolites with 132 patient samples.
This data has been formatted for IntLIM. We load in the data as follows. The bc.ambs.csv meta file contains a list of phenotypic data file, metabolite data file, gene expression data file, metabolite meta file, and gene expression meta file. The function OutputStats will give a summary of the NCI-60 data
inputData <- IntLim::ReadData('bc.ambs.csv',metabid='id',geneid='id')
## Setting options('download.file.method.GEOquery'='auto')
## Setting options('GEOquery.inmemory.gpl'=FALSE)
## [1] "CreateMultiDataSet created"
IntLim::OutputStats(inputData)
## Num_Genes Num_Metabolites Num_Samples_withExpression
## 1 20254 536 108
## Num_Samples_withExpression.1 Num_Samples_inCommon
## 1 132 108
From the OutputStats, we find that we have gene expression data involving 20,254 genes with 108 patient samples and metabolite abundance data involving 536 metabolites with 132 patient samples.
We remove genes with mean belows below the 10th percentile. Furthermore, we remove metabolites with more than 80% missing values. This results in gene expression data involving 18,228 genes and 108 patient samples and metabolite abundance data involving 379 metabolites and 132 patient samples.
inputDatafilt <- IntLim::FilterData(inputData,geneperc=0.10, metabmiss = 0.80)
## [1] "No metabolite filtering by percentile is applied"
IntLim::OutputStats(inputDatafilt)
## Num_Genes Num_Metabolites Num_Samples_withExpression
## 1 18228 379 108
## Num_Samples_withExpression.1 Num_Samples_inCommon
## 1 132 108
We can obtain boxplot distributions of the data as follows. This is used to make figures.
IntLim::PlotDistributions(inputDatafilt)
The principal component analysis is performed on filtered metabolite and gene expression data to obtain visual representations showing how different sub-sections of the data could be grouped into different clusters. Common samples patient samples (either tumor or adjacent non-tumor samples). Blue samples indicate tumor samples and red samples indicate non-tumor samples. Note the clear delineation of samples.
PlotPCA(inputDatafilt, stype = "DIAG", common = T)
## Warning in RColorBrewer::brewer.pal(numcateg, palette): minimal value for n is 3, returning requested palette with 3 different levels
The linear model is for integrating transcriptomic and metabolomics data is: E(m|g,t) = β1 + β2 g + β3 p + β4 (g:p) + ε where ‘m’ and ‘g’ are log-transformed metabolite abundances and gene levels respectively, ‘p’ is phenotype (cancer type, patient diagnosis, treatment group, etc), ‘(g:p)’ is the association between gene expression and phenotype, and ‘ε’ is the error term that is normally distributed. A statistically significant p value of the ‘(g:p)’ association term indicates that the slope relating gene expression and metabolite abundance is different from one phenotype compared to another. We run a linear model on tumor (n = 61) and non-tumor samples (n = 47) that included 18,228 genes and 379 metabolites (total of 6,908,412 possible associations and hence models). For genes and metabolites that had standard deviations of 0 in one of the groups, we assign a p value of NA. The model is run as below by calling RunIntLim(). DistPvalues() allows us to obtain a distribution of p-values for the (g:p) term.
myres <- IntLim::RunIntLim(inputDatafilt, stype="DIAG")
## [1] "Running the analysis on"
##
## NORMAL TUMOR
## 47 61
## [1] "10 % complete"
## [1] "20 % complete"
## [1] "30 % complete"
## [1] "40 % complete"
## [1] "50 % complete"
## [1] "60 % complete"
## [1] "70 % complete"
## [1] "80 % complete"
## [1] "90 % complete"
## user system elapsed
## 347.855 33.232 382.035
IntLim::DistPvalues(myres)
The next step is to process the results of this model by filtering the results of the linear model by FDR-adjusted p-value cutoff (0.10 selected here) for the (g:p) association coefficient and calculate the correlations of the gene-metabolite pairs in each group (tumor and non-tumor) for the filtered results. We further may only interested in results that have an absolute correlation value difference (0.50 selected here). This is done with the ProcessResults function. In addition we also develop a heatmap of the gene-metabolite association correlations for the selected groups.
myres <- IntLim::ProcessResults(myres, inputDatafilt, diffcorr = 0.50, pvalcutoff = 0.05)
IntLim::CorrHeatmap(myres)
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`
dim(myres@corr)
## [1] 2842 4
From this model we find 2842 gene-metabolite correlations that have an association FDR-adjusted p-value of 0.05 and an absolute value correlation difference of 0.5 or greater. The top pairs are shown below.
corr.table <- myres@corr
abs.corrdiff <- abs(myres@corr$NORMAL - myres@corr$TUMOR)
sort.table <- corr.table[order(-abs.corrdiff),]
sort.table[1:10,]
## metab gene NORMAL TUMOR
## 445 adrenate (22:4n6) PDLIM4 0.5925755 -0.6262295
## 1035 docosapentaenoate (n3 DPA; 22:5n3) HIST1H3A -0.5074006 0.6666402
## 454 adrenate (22:4n6) HIST1H3A -0.4935816 0.6429931
## 472 adrenate (22:4n6) GLI3 0.5703712 -0.5602327
## 334 adrenate (22:4n6) GRIA4 0.5279288 -0.6001586
## 1372 guanosine IGFBP4 0.5839669 -0.5426230
## 378 adrenate (22:4n6) IGFBP4 0.5797386 -0.5453728
## 1537 leucine COL14A1 0.5366559 -0.5864622
## 1032 docosapentaenoate (n3 DPA; 22:5n3) PDLIM4 0.5124884 -0.6061952
## 416 adrenate (22:4n6) FHL2 0.5068810 -0.6087784
We can show some example plots of some of these pairs. The first example is the PDLIM4 vs. adrenate (22:4n6). There appears to be an error with the scatterplot.
IntLim::PlotGMPair(inputDatafilt, stype = "DIAG", geneName = "PDLIM4", metabName = "adrenate (22:4n6)")